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Executive  Summary 


The  Virtual  Ship  Architecture  (VSA)  provides  a  framework  for  conducting  distributed 
simulation  in  the  maritime  domain.  It  is  based  upon  the  High  Level  Architecture  (HLA) 
and  supports  the  linking  of  simulations  that  represent  complete  entities,  such  as  warships, 
submarines  and  aircraft,  or  simulations  that  represent  the  systems  of  which  these  entities 
are  composed.  These  include  sensor,  weapon,  countermeasure,  navigation  and  command 
and  control  systems.  It  is  through  linking  these  simulations  that  a  virtual  representation 
of  a  warship  may  be  created. 

An  essential  requirement  for  distributed  simulation  components  to  be  linked  and  operate 
together  is  to  exchange  information  about  the  kinematic  attributes  of  the  entities  that 
are  modelled.  The  kinematic  attributes  include  position,  velocity,  acceleration  and  those 
attributes  that  specify  orientation.  These  attributes  must  be  provided  in  a  consistent 
way,  with  respect  to  a  common  frame  of  reference. 

A  key  part  of  the  Virtual  Ship  Architecture  is  therefore  the  adoption  of  a  common  frame 
of  reference  with  respect  to  which  the  kinematic  attributes  of  entities  are  transmitted 
amongst  the  participating  simulations. 

The  common  frame  of  reference  used  in  the  Virtual  Ship  Architecture  is  the  WGS84 
coordinate  system.  All  entities  in  the  simulation  are  considered  to  have  a  coordinate 
system  fixed  to  them.  The  origin  of  this  coordinate  system  specifies  the  location  of  an 
entity.  The  orientation  of  this  coordinate  system  with  respect  to  the  reference  specifies 
the  entity’s  orientation. 
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1.  Introduction 

The  Virtual  Ship  Architecture  (VSA)  provides  a  framework  for  conducting  distributed 
simulation  in  the  maritime  domain.  It  is  based  upon  the  High  Level  Architecture  (HLA) 
and  supports  the  linking  of  simulations  that  represent  complete  entities,  such  as  warships, 
submarines  and  aircraft,  or  simulations  that  represent  the  systems  of  which  these  entities 
are  composed.  These  include  sensor,  weapon,  countermeasure,  navigation  and  command 
and  control  systems.  It  is  through  linking  these  simulations  that  a  virtual  representation 
of  a  warship  may  be  created.  The  Virtual  Ship  Architecture  Description  Document 
(VSADD)  ill  provides  a  comprehensive  description  of  the  VSA  and  should  be  read  in 
conjunction  with  this  document. 

An  essential  requirement  for  distributed  simulation  components  to  be  linked  and  operate 
together  is  to  exchange  information  about  the  kinematic  attributes  of  the  entities  that 
are  modelled.  The  kinematic  attributes  include  position,  velocity,  acceleration  and  those 
attributes  that  specify  orientation.  These  attributes  must  be  provided  in  a  consistent 
way,  with  respect  to  a  common  frame  of  reference.  Upon  receipt  of  attribute  updates, 
each  federate  will  typically  need  to  transform  these  data  to  coordinate  systems  that  are 
local  to  the  entities  it  represents. 

A  key  part  of  the  Virtual  Ship  Architecture  is  therefore  the  adoption  of  a  common  frame 
of  reference  with  respect  to  which  the  kinematic  attributes  of  entities  are  transmitted 
amongst  the  participating  federates.  This  document  describes  the  common  frame  of 
reference  used  within  the  Virtual  Ship  Architecture  as  a  basis  for  providing  kinematic 
attributes  of  entities. 


2.  The  reference  coordinate  systems 


The  fundamental  reference  coordinate  system  is  a  right-handed  Cartesian  coordinate 
system  whose  origin  is  located  at  the  centre  of  mass  of  the  earth.  The  coordinate  system 
is  considered  fixed  within  the  earth,  so  it  rotates  with  it.  In  this  coordinate  system 
the  z-axis  is  defined  by  the  Conventional  Internation  Origin  (CIO)  which  is  the  mean 
position  of  the  earth’s  rotation  axis  in  the  period  1900-1905  [2].  The  x-axis  is  defined  as 
passing  through  the  mean  Greenwich  meridian  and  the  y-axis  completes  a  right-handed 
coordinate  system.  This  is  illustrated  in  Figure  1. 


This  coordinate  system  is  used  extensively  in  geodesy.  A  principal  use  is  to  define  ellip¬ 
soids  that  approximate  the  shape  of  the  earth.  Typically  these  are  oblate  spheriods  given 
by  the  expression 


2  2  2 

X1  tr  2Z 

- 1_  — - 1 - 

a2  a2  52 


=  1. 


(1) 


In  this  expression  a  is  the  semi-major  axis  and  b  is  the  semi-minor  axis.  Other  quantities 
used  in  reference  to  these  ellipsoids  are  the  flattening 


/  =  ( a-b)/a , 


(2) 


and  the  eccentricity 

e2  =  1  -  b2 /a2  =  2/  -  /2.  (3) 


A  commonly  used  ellipsoid  in  the  military  context  is  that  which  define  the  WGS84  coor- 
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dinate  system  [3].  The  parameters  that  define  this  ellipsoid  are 

a  —  6378137.00000  m, 
b  =  6356752.31425  m, 

1/f  =  298.2572235614, 
e2  =  0.00669437999013. 

In  the  WGS84  coordinate  system  the  position  of  a  point  is  typically  given  as  a  latitude 
4>.  longitude  A  and  height  h  above  the  ellipsoid.  The  meridian  plane  is  defined  as  the 
plane  containing  both  the  z-axis  and  the  point  concerned.  The  normal  to  the  ellipsoid 
passing  through  this  point  lies  in  the  meridian  plane.  The  latitude  is  the  angle  between 
the  normal  to  the  ellipsoid  and  the  x  -  y  plane,  measured  in  the  meridian  plane.  The 
longitude  is  the  angle  between  the  Greenwich  meridian  and  the  meridian  plane  of  the 
point.  The  height  coordinate  is  the  height  of  the  point  above,  or  below,  the  ellipsoid, 
measured  along  the  normal.  Typically  the  latitude  is  given  as  <j>  G  [— 7t/2,  7t/2]  and  the 
longitude  is  given  as  A  €  (— 7T,  7t]  .  Latitude  is  measured  positive  to  the  North  of  the 
equator  and  negative  to  the  South.  The  longitude  is  measured  negative  to  the  West  and 
positive  to  the  East  of  the  Greenwich  meridian. 


f’ 


Figure  1:  The  G  coordinate  system. 

In  the  remainder  of  this  document  this  coordinate  system  is  central.  We  denote  it  as  the 
G  coordinate  system.  In  addition,  we  will  use  superscripts  to  denote  quantities  that  are 
given  with  respect  to  this  coordinate  system.  For  example,  the  position  vector  of  a  point 
given  with  respect  to  this  coordinate  system  will  be  donted  as  rG. 
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For  a  point  positioned  on,  or  near,  the  Earth’s  surface  it  is  often  appropriate  to  work 
with  a  local  coordinate  system.  Let  us  define  a  coordinate  system  about  some  point.  We 
will  call  this  the  local  geodetic  coordinate  system,  and  denote  it  as  LG.  It  is  illustrated 
in  Figure  2.  The  zLG-axis  is  defined  by  the  exterior  normal  to  the  ellipsoid  that  passes 
through  the  point,  the  yLG~axis  is  perpendicular  to  this  and  defines,  with  the  zLG-axis, 
a  plane  that  contains  the  zG  axis.  The  yLG- axis  is  said  to  define  geodetic  North.  The 
xLG-axis  is  chosen  to  complete  a  right  handed  coordinate  system  and  define  East.  The 
transformation  between  the  G  and  LG  coordinate  system  is  given  in  Appendix  B.  Note 
that  the  LG  coordinate  system  may  be  defined  at  any  point.  It  need  not  have  its  origin 
on  the  surface  of  the  ellipsoid. 


Figure  2:  The  LG  coordinate  system. 

3.  Specifying  the  kinematic  attributes  of  an  entity 

Within  the  Virtual  Ship  Architecture  are  notions  of  entities  that  occupy  a  finite  volume 
of  space  and  these  include  the  CompositeEntity,  ComponentEntity  and  Sensor  Task.  An 
instance  of  a  CompositeEntity  object  represents  a  complete  entity  within  the  VSA.  Ships, 
submarines,  aircraft,  missiles  and  countermeasures  are  all  examples  of  CompositeEntity 
objects.  An  instance  of  a  ComponentEntity  object  represents  entities  that  are  compo¬ 
nents  of  a  CompositeEntity.  Sensor,  weapon,  countermeasure,  navigation  and  command 
and  control  systems  are  all  examples  of  ComponentEntity  objects.  An  instance  of  the 
SensorTask  object  represents  a  particular  task  being  performed  by  a  sensor,  such  as  a 
radar,  ESM  or  sonar  and  is  characterised,  in  part,  by  the  sensor  beam. 
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To  specify  the  kinematics  of  each  of  these  types  of  entity,  it  is  assumed  that  they  have 
associated  with  them  a  right-handed  Cartesian  coordinate  system.  The  x-axis  of  this 
coordinate  system  is  typically  in  the  direction  that  would  commonly  be  identified  as 
forwards,  sometimes  known  as  the  figure  axis.  The  z-axis  is  in  the  direction  typically 
identified  as  up  and  the  y-axis  completes  a  right-handed  coordinate  system.  The  origin 
of  the  coordinate  systems  associated  with  CompositeEntity  and  ComponentEntity  objects 
is  located  at  the  centre  of  mass  of  the  entity. 

For  notational  purposes,  the  coordinate  system  associated  with  a  CompositeEntity  object 
instance  shall  be  denoted  by  E,  that  associated  with  a  ComponentEntity  object  instance 
shall  be  denoted  by  e  and  that  associated  with  a  sensor  task  shall  be  denoted  by  S.  These 
coordinate  systems  are  illustrated  in  Figures  3  and  4. 


Figure  3:  The  E  and  e  coordinate  systems. 

To  specify  the  kinematics  of  these  entities  requires  that  the  motion  of  the  origin  of  the 
coordinate  system  be  given,  and  the  orientation  of  the  coordinate  system  given  with 
respect  to  some  reference  frame.  The  motion  of  the  origin  is  given  as  a  position,  velocity 
and  acceleration  vector.  The  notations  o% \ ,  o%\  and  b%\  are  used  to  specify  these  values. 
The  subscript  Cl  denotes  the  coordinate  system  of  which  this  point  is  the  origin  and 
the  superscript  02  denotes  the  coordinate  system  within  which  the  origin  is  specified. 
The  orientation  is  specified  through  the  use  of  Euler  angles  that  define  the  rotation 
transformation  from  the  reference  coordinate  system  C%  to  the  entity  coordinate  system 
C\ .  The  convention  for  defining  the  Euler  angles  is  described  in  Appendix  C  and  the  set 
of  three  Euler  angles  is  denoted  by  QC2_>Cl .  In  addition  to  specifying  the  Euler  angles 
their  time  derivatives  may  also  be  given  to  specify  the  rotation  rate  of  the  entity,  and 
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these  are  denoted  as  tiC2~"Cl .  The  relationship  between  time  derivatives  of  the  Euler 
angles  and  the  angular  velocity  of  a  rigid  body  is  given  in  Appendix  F. 


4.  Object  attributes  that  convey  coordinate  data 

Within  the  VS-FOM  there  are  a  number  of  object  attributes  that  convey  kinematic 
information.  The  remainder  of  this  section  describes  these  attributes. 

4.1  The  CompositeEntity  object  class 

The  composite  entity  has,  among  others,  the  following  attributes: 

Position, 

Velocity, 

Acceleration, 

Orientation, 

OrientationRate. 

The  Position  attribute  is  the  position  vector  of  the  origin  of  the  E  coordinate  system, 
given  with  respect  to  the  G  system.  This  is  written  as 

Position  =  o°. 
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The  Velocity  attribute  is  the  velocity  of  the  origin  of  the  E  coordinate  system,  given  with 
respect  to  the  G  system.  This  is  written  as 

Velocity  =  of. 

The  Acceleration  attribute  is  the  acceleration  of  the  origin  of  the  E  coordinate  system, 
given  with  respect  to  the  G  system.  This  is  written  as 

Acceleration  =  of. 

The  Orientation  attribute  is  the  Euler  angle  set  that  defines  the  rotation  transformation 
from  the  G  system  to  the  E  system.  This  is  written  as 

Orientation  =  QG~*E. 

The  OrientationRate  attribute  is  the  time  derivative  of  the  Euler  angle  set  that  defines 
the  rotation  transformation  from  the  G  system  to  the  E  system.  This  is  written  as 

OrientationRate  =  0G_*E. 


4.2  The  ComponentEntity  object  class 

The  component  entity  has,  among  others,  the  following  attributes: 

RelativePosition, 

RelativeOrientation. 

The  RelativePosition  attribute  is  the  position  vector  of  the  origin  of  the  e  coordinate 
system,  given  with  respect  to  the  E  system  of  the  CompositeEntity  object  of  which  this 
ComponentEntity  is  a  part,  otherwise  known  as  the  parent  entity.  This  is  written  as 

RelativePosition  =  of. 

The  RelativeOrientation  attribute  is  the  Euler  angle  set  that  defines  the  rotation  trans¬ 
formation  from  the  E  system  to  the  e  system.  This  is  written  as 

RelativeOrientation  =  QE~*e. 


4.3  Sensor  tasks 

The  SensorTask  object  instance  has,  among  others,  the  following  attributes: 

BeamPatternReference, 

BeamPatternOrientation, 

BeamPatternOrientationRate. 

The  origin  of  the  beam  pattern  coordinate  system  S  is  assumed  to  be  colocated  with 
the  origin  of  the  coordinate  system  attached  to  the  ComponentEntity  with  which  the 
SensorTask  is  associated.  The  BeamPatternReference  attribute  specifies  the  coordinate 
system  with  respect  to  which  the  orientation  of  the  S  coordinate  system  is  given.  The 
BeamPatternOrientation  attribute  is  the  Euler  angle  set  that  defines  the  rotation  trans¬ 
formation  from  the  reference  system  to  the  S  system.  The  BeamPatternOrientationRate 
attribute  provides  the  time  derivative  of  these  Euler  angles. 
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If  BeamPatternReference 

and 


=  1  then  the  reference  system  is  the  G  system,  so 
BeamPatternOrientation  =  QG~*S, 


BeamPatternOrientationRate  =  QG'S. 


If  BeamPatternReference  =  2  then  the  reference  system  is  the  LG  system,  so 

BeamPatternOrientation  =  QLG^S, 


and 


BeamPatternOrientationRate  =  QLG^S. 


If  BeamPatternReference  =  3  then  the  reference  system  is  the  E  system,  that  is  the 
coordinate  system  of  the  CompositeEntity  object  instance  that  is  the  parent  of  the  Com- 
ponentEntity  with  which  the  task  is  associated.  Hence 

BeamPatternOrientation  =  GE^S, 


and 


BeamPatternOrientationRate  =  GE  s 


The  LG  system  will  be  appropriately  used  as  the  reference  in  the  case  of  a  stabilised 
sensor  and  the  E  system  might  be  used  otherwise. 

4.4  Tracks 

The  Track  class  and  its  subclasses  have  various  attributes  that  convey  kinematic  infor¬ 
mation.  A  track  is  considered  to  represent  a  point  in  space  and  does  not,  therefore, 
require  specification  of  orientation.  In  the  first  instance,  consider  the  attributes  of  the 
AbsoluteTrack  subclass.  Those  that  convey  kinematic  information  are 

Position, 

PositionError, 

Velocity, 

VelocityError. 

The  Position  attribute  is  the  position  vector  of  the  track  given  with  respect  to  the  G 
system.  The  Velocity  attribute  is  the  velocity  of  the  track  given  with  respect  to  the  G 
system.  The  error  attributes  are  measured  with  respect  to  the  same  coordinate  system 
as  the  quantities  to  which  they  refer. 

The  attributes  of  the  RelativeTrack  subclass  that  convey  kinematic  information  are 

FootpointPosition, 

FootpointPositionError, 

Bearing, 

BearingError, 

Elevation, 

ElevationError, 

Range, 

RangeError, 

RangeRate, 

RangeRateError. 
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The  FootointPosition  attribute  is  the  position  vector  of  the  point  with  respect  to  which 
the  track  is  specified,  given  with  respect  to  the  G  system. 

The  Bearing  attribute  is  the  bearing  to  the  track,  given  with  respect  to  the  LG  coordinate 
system  defined  at  the  footpoint.  The  bearing  is  the  angle  between  the  xLG-axis  and  the 
vertical  plane  that  contains  the  zLG-axis  and  the  track  position.  This  is  illustrated  in 
Figure  5. 


ZLG  (Vertical) 


(North) 


Figure  5:  The  bearing  and  elevation  attributes  of  the  RelativeTrack  object  class. 

The  Elevation  attribute  is  the  elevation  to  the  track  point,  given  with  respect  to  the  LG 
coordinate  system  defined  at  the  footpoint.  The  elevation  is  the  angle  between  a  line 
joining  the  footpoint  to  the  track  and  the  xhG  —  yLC  plane. 

The  Range  attribute  is  the  straight-line  distance  between  the  footpoint  and  the  track. 
The  RangeRate  attribute  is  the  time  derivative  of  this  quantity. 

The  error  attributes  are  measured  with  respect  to  the  same  coordinate  system  as  the 
quantities  to  which  they  refer. 


5.  Dead- reckoning 

To  reduce  the  frequency  of  certain  kinematic  attribute  updates  the  technique  of  dead¬ 
reckoning  is  used.  Dead-reckoning  involves  extrapolating  the  kinematic  attributes  of  an 
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entity  in  time.  A  federate  responsible  for  updating  the  kinematic  attributes  of  an  entity 
will  typically  also  dead-reckon  them  and  provide  attribute  updates  when  the  deviation 
between  the  actual  and  dead-reckoned  values  exceeds  some  desired  accuracy. 

To  elaborate  by  way  of  an  example,  suppose  a  federate  provides  the  following  at  time  t\: 

ol{h),  of(ti),  oG(t 1), 

that  is  the  position,  velocity  and  acceleration  of  a  composite  entity.  Extrapolating,  the 
dead-reckoned  position  at  time  t  is 

DRO°(t)  =  o°(ti)  +  o°(<i)(t  -  h)  +  50°(ti)(i  -  h)2.  (4) 

If  the  position  is  required  to  be  known  throughout  the  federation  with  accuracy  e,  then 
an  attribute  update  is  required  only  when 

|o°(t)  -dr  of  (t)  |  >  e.  (5) 

The  significant  advantage  of  this  approach  is  that  if  the  acceleration  of  an  entity  varies 
slowly  then  a  large  time  will  elapse  before  (5)  is  satisfied  and  an  attribute  update  is 
required. 

There  are  a  number  of  approaches  to  dead-reckoning.  These  primarily  concern  which  of 
the  attributes  should  be  extrapolated  in  time  and  whether  they  should  be  extrapolated 
to  first  or  second  order  in  time.  The  CompositeEntity  attribute  DRAlgorithm  provides 
an  enumeration  of  possible  algorithms.  For  each  object  instance,  the  federate  updating 
the  kinematic  attributes  of  the  object  will  determine  which  algorithm  is  appropriate. 
Federates  that  subscribe  to  these  attributes  will  require  the  capability  to  utilise  any  of 
the  dead-reckoning  algorithms. 

The  enumerators  One,  Two,  Three,  Four,  Five  and  Six  specify  which  dead-reckoning 
algorithm  to  use  and  these  are  detailed  in  the  following.  It  is  assumed  that  the  most 
recent  attribute  update  has  occurred  at  time  t\. 

DRAlgorithm  =  One 

<%(t) 
o°{t) 

*i(t) 
nc^s(t) 

QG~>E(t) 

DRAlgorithm  =  Two 

o°E(t ) 

6g(t) 

ot(t) 

QG^E(f) 

Cic-E(t) 


=  <>e(0  +  of(£i)(i-ti), 


aG 


(*l)> 


DG^E(ti), 


(7) 


=  oG(ii), 

=  oG(fx) 

=  oG(ti),  (6) 

=  ^E(t1), 

=  fiG-E(ti). 
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DRAlgorithm  =  Three 

o£(*)  =  °!0i)  +  -  ti)  +  -  ti)2, 

Oe0)  =  bG(t i)  +  o£(*i)(i  -  ti), 

«eW  =  ag(fi), 

ftG_+E(f)  = 

ciG^E(t)  =  f2G->E(ti). 


DRAlgorithm  =  Four 


O  eW  =0G(ti), 
o£(*)  =  o£(*i), 

OeW  =  o|(ti), 

QG_>E(t)  =  nG^E{h)  +  nG^E(ti)(t  -  ti), 

nG~*E(t)  =  CiG~"E{ti). 


DRAlgorithm  =  Five 


o£(i)  =  oG(ti)  +6G(ti)(t-ti), 

OeW  =  Oe(<i)» 

6gW  =  5l(*i), 

QG^E  (t)  =  DG_>E(ti)  +  flG_*E(fi)(f  —  ti), 

QG~*E(f)  =  flG^E(ti). 


(8) 


(9) 


(10) 


DRAlgorithm  =  Six 

°eW  =  °e(^i)  +  —  h)  +  50g(ti)(t  —  ti)2, 

o°(t)  =  6|(ti)  +  oG(ti)(t  -  ti), 

«£(*)  =  °£(*i),  (11) 

DG-+E(f)  =  fiG^E(ti)  +  DG^E(ti)(t  -  ii), 

DG“*E(t)  =  DG^E(ti). 


Dead-reckoning  is  also  used  to  extrapolate  the  orientation  of  the  coordinate  system  asso¬ 
ciated  with  the  beam  pattern  that  defines  the  SensorTask  object.  The  attribute  Beam- 
PatternDRAlgorithm  provides  an  enumerated  description  of  the  algorithm  used.  In  the 
following  it  is  assumed  that  the  coordinate  system  S  is  given  with  respect  to  the  G  system. 
The  formulae  are  the  same  in  the  event  that  the  LG  or  E  systems  are  used  as  reference. 

BeamPatternDRAlgorithm  =  One 

=  nG-s(h): 

nG-s(t)  =  ciG~*s(ti). 


BeamPatternDRAlgorithm  =  Two 

fiG_>s(f)  =  nG^s(ii)  +  DG^s(fi)(f  —  ti), 
f2G_*s(f)  =  £lG~*s(ti). 


(13) 
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Appendix  A  -  Transforming  between 
Cartesian  coordinates  and 
latitude,  longitude  and  height 

The  transformation  from  a  representation  of  a  point  as  a  latitude,  longitude  and  height 
(0,  A ,h)  to  Cartesian  coordinates  (x,y,z)  in  the  G  system  is  straight  forward.  We  have 

x  =  (N  +  h)  cos  0  cos  A.  (Al) 

y  =  (N  +  h)  cos  0  sin  A,  (A2) 

z  —  ((1  -  e2)Ar  +  h)  sin0,  (A3) 

where  _ 

N  =  a/yj  1  —  e2  sin2  0.  (A4) 


The  inverse  transformation  is  a  little  more  complex.  From  (Al)  and  (A2)  define 

p  =  y/{x2  +  y2)  =  (N  +  h)  cos  0. 
noting  that  —  7r/2  <  0  <  0/2.  Using  this  in  (A3)  yields 

2  =  (p/  cos  0  —  e2 N)  sin  0. 


If  we  define 


f  (0)  =  p  tan  0  —  e2N  sin  0  —  2, 

then  0  is  a  root  of  /  and  may  be  found  using  appropriate  root  finding  techniques. 


(A5) 

(A6) 

(A7) 


The  Newton-Raphson  method  is  a  technique  that  works  well  in  this  case.  To  begin  the 
iterations  using  this  method  requires  an  initial  estimate  of  the  root  which  can  be  obtained 
as  follows.  Using  (A5)  in  (A3)  yields 

tan0  ^1  -  =  ztv -  (A8) 

Now  N  is 'of  the  order  of  the  radius  of  the  Earth,  so  for  objects  for  which  h  <S  AT, 
N/(N  +  h)  «  1  and  (A8)  gives  a  first  approximation  to  the  latitude  as 


0^°^  =  arctan 


p(l  -  e2) 


(A9) 


It  should  be  noted  that  this  approximation  is  of  high  accuracy  and  only  one  or  two 
interations  of  the  Newton-Raphson  method  are  typically  required. 


Having  determined  0  it  remains  to  determine  A  and  h.  A  convenient  way  to  get  A  is  to 
first  note  that 


A  sin(A/2)  _  2sin(A/2)  cos(A/2)  _  sin  A 
tan  2  =  cos(A/2)  ~  2cos2(A/2)  _  cosA+  1' 


(A10) 


Using  (Al),  (A2)  and  the  definition  of  p  yields 


y 

x  +  p' 


(All) 


so 


A  =  2  arctan 

h  is  obtained  directly  from  (A5)  as 


y 


x+p 
h  =  p/  cos  0  —  N. 


(A12) 

(A13) 
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Appendix  B  -  Transforming  between  the 
G  and  LG  coordinate  systems 


The  unit  vectors  that  define  the  LG  coordinate  system  are  given  in  terms  of  the  latitude 
and  longitude  of  the  origin  of  this  system.  Let  this  point  have  Cartesian  coordinates 
(. xG,yG,zG )  with  respect  to  the  G  coordinate  system,  which  may  be  alternatively  given 
as  (0,  A,  h).  With  reference  to  Figure  2,  we  have 

iLG  =  —  sinAfG  +  cosAjG, 

jhG  =  —  sin0cosAiG  —  sin  0  sin  XjG  +  cos  0fcG,  (Bl) 

kLG  =  cos  0  cos  A  iG  +  cos  0  sin  A  jG  +  sin  0  kG . 


If  eLG  denotes  a  unit  vector  with  components  given  with  respect  to  the  LG  coordinate 
system,  and  eG  denotes  the  same  vector  with  components  given  with  respect  to  the  G 
system  then 

eG  =  i?LG_>GeLG,  (B2) 

where  the  rotation  matrix  -RLG— >G  is  given  as 


R 


LG- 


—  sin  0  cos  A 

—  sin  0  sin  A 

COS0 


cos  0  cos  A ' 
cos  0  sin  A 
sin  0 


(B3) 


From  (B2)  we  also  have 


—  i?-1  eG  =  Rt  ec 

~  -rlLG— <-Ge  — >G  c 


-Rg— ►  LGe<" 


(B4) 


where  the  property  of  rotation  matrices  that  their  inverse  is  equal  to  the  transpose  gives 
us 


or  writing -it  out  in  full 


Rc 


Rg^  LG  — 

-“LG— >G> 

(B5) 

—  sin  A 

cos  A 

°  \ 

—  sin  0  cos  A 

—  sin  0  sin  A 

COS  0 

(B6) 

cos  0  cos  A 

cos  0  sin  A 

sin0  / 

Consider  now  the  transformation  of  position  vectors  give  with  respect  to  the  G  and  LG 
coordinate  systems.  Denote  by  oGG  the  position  vector  of  the  origin  of  the  LG  system, 
given  with  respect  to  the  G  coordinate  system.  It  is  assumed  that  the  origin  of  the  G 
system  is  the  null  vector.  Let  rG  denote  the  position  vector  of  a  point  with  respect  to 
the  G  coordinate  system  and  let  rLG  denote  the  position  vector  of  the  same  point  with 
respect  to  the  LG  system.  Then  we  have 

+  (B7) 

and 

r“  =  -  d»,).  (Bg) 

where  the  relationship  between  f?G_> LG  and  Rlg-+g  has  been  noted  in  (B5). 

It  is  important  to  bring  out  the  fact  that  the  elements  of  these  rotation  matrices  are 
dependent  upon  oGG,  whether  the  position  of  this  point  is  specified  as  a  latitude,  lon¬ 
gitude  and  height,  or  in  Cartesian  coordinates.  It  is  noted  that  the  rotation  matrices 
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are  expressed  more  cleanly  in  terms  of  (<p,  A ,h),  but  these  may  be  obtained  in  terms  of 
(xG,yG,zG)  through  the  transformation  defined  in  Appendix  A. 


When  transforming  the  velocities  between  the  G  and  LG  systems  the  time  derivative  of 
the  matrices  Ag-.lg  and  RLg->g  are  required.  The  time  derivative  of  Rg->lg  is 


Rg 


LG  —  (f) 


0 

—  cos  4>  cos  A 

—  sin  d>  cos  A 


(—  cos  A 
sin  6  sin  A 
—  cos  &  sin  A 


0 

—  cos  <f>  sin  A 

—  sin  sin  A 
—  sin  A 

—  sin  <f>  cos  A 
cos  <±>  cos  A 


°  \ 

—  sin  <f>  ] 

COS  4>  ) 


(B9) 


The  time  derivative  of  RLg-^g  is  the  transpose  of  this. 
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Appendix  C  -  The  rotation  matrix  in 
terms  of  the  Euler  angles 

The  rotation  transformation  from  one  frame  to  another  is  considered  as  composed  of 
three  rotations,  denoted  by  fl  =  (flx,  fly,  ez).  The  convention  for  defining  these  rotations 
and  the  expression  for  the  rotation  matrix  are  defined  in  this  appendix. 


Figure  Cl:  The  Euler  angles. 


Let  r1  denote  the  components  of  a  vector  with  respect  to  the  coordinate  system  denoted 
by  I,  and  suppose  that  after  transformation  through  a  rotation  described  by  the  Euler 
angle  set  fl  the  components  of  this  vector  are  rR.  The  first  rotation  is  flz  about  the  zr 
axis,  measured  positive  as  clockwise  when  looking  along  the  axis  and  illustrated  in  Figure 
Cl.  All  other  rotations  are  defined  using  this  convention.  The  matrix  of  this  rotation  is 

(cos  flz  sin  f2z  0  \ 

—  sin  flz  cos  flz  0 

0  0  1/ 


The  second  rotation  is  fly  about  the  y(-axis,  as  illustrated  in  Figure  Cl.  Here  the  super¬ 
script  '  indicates  the  coordinate  system  obtained  after  the  first  rotation.  The  matrix  of 


this  rotation  is 


COS  Chy  0  sin  fly 


sin  fly  0  cos  fly 


The  final  rotation  is  flx  about  the  yf-axis,  as  illustrated  in  Figure  Cl.  The  superscript 
"  indicates  the  coordinate  system  obtained  after  two  rotations.  The  matrix  of  this  final 
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rotation  is 

(l  0 
0  cos  flx 
V  0  —  sin  Qx 

Composing  the  matrices  gives  the  transformation  as 


rR  =  0 


0 

cos 

—  sin  fL 


Evaluating  the  matrix  product  gives  the  rotation  matrix  as 


sin  Qz 
cosfiz 
0 


cos  Qy  cos  f 2- 
—  cos  Qx  sin  Q~  +  sin  Gx  sin  Gy  cos  G2 
,  sin  Gz  sin  G;  +  cos  Gx  sin  Gy  cos  Gz 


R  = 


(Cl) 


cos  sin  G~  —  sin  Gy 

cos  Gx  cos  Gc  +  sin  Gx  sin  Gy  sin  Gz  sin  Gx  cos  Gy 
—  sin  Gx  cos  G-  +  cos  flx  sin  Gy  sin  Gz  cos  Gz  cos  Gy , 


The  time  derivative  of  the  rotation  matrix  is 


u  u 

sin  Gx  sin  G2  +  cos  Gx  sin  Gy  cos  G2  —  sin  Gx  cos  G2  +  cos  Gx  sin ! 
cos  Gx  sin  G;  —  sin  Gx  sin  Gy  cos  G2  —  cos  Gx  cos  G-  —  sin  Gx  sin ! 

-  sin  Gy  cos  Gz  —  sin  Gy  sin  Gz 

Gx  cos  Gy  cos  Gx  sin  Gx  cos  Gy  sin  Gz 
;  Gx  cos  Gy  cos  G2  cos  Gx  cos  Gy  sin  Gz 

—  cos  Gy  sin  Gz  cos  Gy  cos  G; 

—  cos  Gx  cos  G-  —  sin  Gz  sin  Gy  sin  Gz  —  cos  Gx  sin  Gz  +  sin  Gx  sin  Gy  cos  Gz 
sin  Gx  cos  Gx  —  cos  Gx  sin  Gy  sin  Gz  sin  Gx  sin  Gz  4-  cos  Gx  sin  Gy  cos  Gz 


cos  Gx  sin 
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Appendix  D  -  Recovering  the  Euler  angles 
from  a  rotation  matrix 


Suppose  we  have  the  rotation  matrix  R  with  elements 


/  Rn  R\2  R\z  \ 

R  =  Rn  R22  R23  ■ 

(Dl) 

Equating  (Dl)  and  (C2) 

\  Rzi  R32  R3Z  J 

we  have 

sin  fly  —  —Ri3- 

(D2) 

We  may  choose  Qy  G  [— 7r/2,7r/2]  so 

cos  Qy  =  (1  —  i?i3)1//2. 

(D3) 

Provided  |i?i3|  7^  1,  that 

is  Q y  7^  ±7t/2,  we  have 

COS  =  i?n/(l  -  ^13)1/2) 

(D4) 

sinfi2  =  iW(l-J??3)1/2, 

(D5) 

sin  flx  —  #23/(1  -  -fti3)1/2> 

(D6) 

cos  flx  =  Rzz/{  1  -  R2iz)1/2- 

(D7) 

If  we  define  the  arccos  function  on  the  interval  [0, 7r]  then  the  sign  of  the  sin  function 
indicates  whether  the  angle  is  positive  or  negative.  From  (Dl)  and  (C2)  we  deduce  that 
sgn(sinffy)  =  sgn(i?i2)  and  sgn(sinffy)  =  sgn(.R23).  Thus 

Dz  =  sgn(I?23)  arccos  (i?33/ (1  -  Rj3)1/7),  (D8) 

Clz  =  sgn(i?i2)  arccos(i?n / (1  -^i3)1/2)-  (D9) 

For  completeness,  consider  the  case  where  |i?i3|  =  1,  so  Qy  =  ±tt/2.  In  the  event  that 
Qy  =  7r / 2  the  rotation  matrix  is 

/  0  0  — 1  \ 

R  =  I  sin(ffy  —  fl2)  cos(Qx  —  ffy)  0  I  .  (DIO) 

Vcos(fIz  -  D2)  -  sin(ffy  -  ffy)  0  ) 

From  this  we  deduce  that 

ffy  —  Qz  =  sgn(I?2i)  arccosi?22-  (Dll) 

Hence  ffy  and  Qz  can  be  chosen  to  be  any  values  that  satisfy  (Dll).  This  follows  from 
the  fact  that  for  a  coordinate  system  that  is  derived  from  another  through  a  sequence  of 
rotations  that  includes  f \y  =  7r/2,  there  is  no  unique  way  of  expressing  the  transformation 
in  terms  of  the  Euler  angles.  However,  in  practical  circumstances  there  may  be  other 
criteria  that  determine  the  values  chosen  for  Slx  and  Qz.  Continuity  of  the  Euler  angles 
may  be  significant  in  this  regard. 

In  the  case  £ly  =  —  -jt/2  the  rotation  matrix  becomes 

/  o  on 

R=  I  -sin(ffy  +  D2)  cos(ffy  +  fl2)  0  j,  (D12) 

\-cos(Qz  +  D2)  -sin(ffy  +  ffy)  0/ 

from  which  we  deduce 

=  -sgn(i?i2)  arccos  R2 2-  (D13) 
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Appendix  E  -  Recovering  the  time  derivative 


of  the  Euler  angles  from  the  time 
derivative  of  a  rotation  matrix 

Taking  the  time  derivative  of  (Dl)  and  equating  with  (C3)  yields 

Rn  =  —Cly  sin  Qy  cos  Qz  —  Qz  cos  Qy  sin  (El) 

R12  =  —Qy  sin  Qy  sin  Qz  +  Q-  cos  Qy  cos  Qz,  (E2) 

Ris  — -QyCosQy.  (E3) 

i?23  =  Qx  cos  Qx  cos  Qy  —  Qy  sin  Qx  sin  Qy,  (E4) 

/Z33  =  —  sin  cos  Qy  —  Qy  cos  Qx  sin  Qy.  (E5) 

Assume  that  cosf2y  ^  0,  that  is  Qy  =  ±7t/2.  From  (E3) 

Qy  —  -R\z/ cosQy.  (E6) 

Multiplying  (El)  by  sinfl2  and  (E2)  by  cosQ2  and  taking  the  difference  yields 

Qz  =  (i?i2  cos  Qz  —  R11  sin  Qz)  j  cos  Qy.  (E7) 

Similarly,  multiplying  (E4)  by  cosflx,  (E5)  by  sinfl^  and  taking  the  difference  yields 

Qx  =  (JR23  cos  Qx  —  H33  sin  Qx)/  cos  Qy .  (E8) 


Assembling  these  results,  provided  cos  Qy  ^  0, 

fir  =  (i?23  cos  Qx  —  R33  sin  Qx )/  cos  Qy, 

Qy  —  —R\z/  cosQy,  (E9) 

=  (i?i2  cos Qz  —  R\\  sm.Qz)/ cosfl^. 

For  completeness,  consider  the  case  where  cosfly  =  0,  so  Qy  =  ±7t/2.  (E1)-(E5)  become 

R\  1  =  —Qy  sin  Qy  cos Qz, 

R 12  =  —Qy  sin  Qy  sin  Qz, 

Riz  —  0,  (E10) 

R23  —  —Qy  sin  Qx  sin  Qy, 

Rz3  =  —  Qy  cosflxSinQ^. 

From  these  we  can  deduce  the  following  equivalent  expressions  for  Qy: 


Qy  =  —(R n  cosQz  +  Ri2  sinfiz)/ sinff^, 

(Ell) 

Qy  =  -(R.23  sin  Qx  +  Rzz  cosff^)/ sinQ^, 

(E12) 

noting  that  sinfly 

=  ±1. 
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To  obtain  expressions  for  Qx  and  Qz  we  need  to  consider  the  expressions  for  the  time 
derivatives  of  the  remaining  elements  of  the  rotation  matrix.  We  have  in  the  case  S2y  = 
7t/2 

R2I  ~~  (S2X  S22)  COS (S2X  Sl2), 


R-22  —  (S2X  )  Sm(Qx  S"l2), 
R31  —  (S2X  S72)  sin(S"2x  Sl2), 


(E13) 


R32  —  (S2*  Sl2)  cos(S2x  S22). 


In  the  case  Qy  =  —  ir/2 

R21  —  —  (^x  +  Si2)  cos(S2x  +  S22), 
R22  =  -(Six  +  S2*)  sin(Qx  +  Sl2), 
jF?3i  =  (^x  +  Si*)  sin(Qx  +  Sl2), 
R32  =  —(Six  +  Si*)  cos(S"2x  +  S22). 


(E14) 


Note  that  we  can  substitute  for  sin(Qx  —  Sl2)  and  cos(S2x  —  Clz)  using  (DIO)  and  for 
sin(fix  +  fi2)  and  cos(Dx  +  S22)  using  (D12).  Hence  we  get  for  the  case  Qy  =  7t/2 


R21  —  ^(six  -  si*), 

R22  =  —  i?2l(six  —  sir), 

R31  =  ^32  (six  —  si*), 
R32  =  — -R3i(six  —  si*), 

and  in  the  case  fly  =  —  n/2 

R21  =  — i?22(six  +  si*), 
R22  =  #21  (six  +  six), 
R31  —  — -R32(six  +  si*), 
R32  =  i?3i(six  +  si*). 


(E15) 


(E16) 


Given  that  R  is  a  rotation  matrix,  not  all  the  coefficients  on  the  right  hand  side  of  (El 5) 
and  (El  6)  can  be  zero  so  a  unique  solution  can  be  obtained  for  Six  —  Si*  in  the  case  where 
Q.y  =  7r/2  and  for  Six  +  Si2  in  the  case  where  f ly  =  -tt/2.  As  for  the  case  of  determining 
the  Euler  angles  given  the  rotation  matrix,  the  freedom  in  solution  for  f2x  and  S~2*  may 
be  exploited  and  perhaps  used  to  ensure  continuity  of  S2X  and  S72. 
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Appendix  F  -  The  relationship  between  the 
time  derivative  of  the  Euler  angles 
and  the  angular  velocity 

Consider  two  Cartesian  coordinate  systems.  The  first  is  stationary  and  the  second  is 
in  rotation  with  respect  to  it.  At  any  instant,  the  transformation  from  the  first  frame 
to  the  second  is  described  by  a  rotation  matrix  defined  in  terms  of  the  Euler  angles  in 
the  convention  used  previously.  The  Euler  angles  are  functions  of  time.  The  rotation 
of  the  second  coordinate  system  may  be  equivalently  described  by  the  angular  velocity 
of  this  frame  relative  to  the  first.  Let  vl  denote  the  velocity  of  a  point  with  respect  to 
the  first  frame.  Let  uR  denote  the  velocity  of  this  same  point  with  respect  to  a  second 
coordinate  system  that  rotates  with  the  angular  velocity  of  the  second  coordinate  system 
define  above,  but  which  is  instantaneously  cooincident  with  the  first  frame.  We  have 

vl  —  vR  +  u  x  r1.  (Fl) 

where  u>  is  the  angular  velocity  of  the  rotating  coordinate  system  measured  with  respect 
to  the  first  coordinate  system,  r 1  is  the  position  vector  of  the  point  given  with  respect  to 
the  first  coordinate  system,  which  is  instantaneously  the  same  as  it  would  be  given  with 
respect  to  the  rotating  system. 


Let  rR  be  the  position  vector  of  the  same  point  given  with  respect  to  the  second  coordinate 
system.  We  have 

rR  =  (F2) 

where  is  the  rotation  matrix  that  associates  the  coordinate  systems.  Differentiating 
with  respect  to  time  we  have 

rR  =  R^r1  +  R^Rr\  (F3) 

In  the  second  frame  of  reference  the  point  is  stationary  so  rR  =  0  and 

f1  =  -R^R^r1  =  u  x  r1.  (F4) 


Now 

u  x  r1  =  ( uyz 1  -  uzyl)  i  +  ( uzxl  -  uxzl)  j  4-  {u)xyl  -  uyxl)  k , 
which  can  be  written  in  matrix  form  as 


Using  this  is  (F4)  gives 


Carrying  out  the  differentiation  we  get 


(  0 

UZ 

~uv\  .  / 

f  0 

1 

°\ 

-uz 

0 

LOx  1  =  flz 

-1 

0 

0 

V  Uy 

~0JX 

0  J  ' 

V  0 

0 

0  J 

+  fly 

+ 


/  0  0  -  cos  flz 

10  0  —  sin  2 

\  cos  Qz  sin  Qz  0 


(  0 

sin  fly 

\  cos  fly  sin  flz 


—  sin  fly 
0 

—  COS  fly  cos  fl~ 


-  cos  fly  sin  Q 
cos  fly  sin  flz 


(F5) 

(F6) 


(FT) 


(F8) 
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Hence 


L<0X  —  fix  COS  fly  COS  fl^  fly  fl z  5 

ujy  =  fix  cos  fly  sin  Qz  +  fly  cos  flz,  (F9) 

Ct^2  —  fl,,  fix  sin  fly . 

This  recovers  the  formulae  given  in  [4]  (page  610). 

It  is  often  required  to  have  the  components  of  the  angular  velocity  with  respect  the  the 
rotating  coordinate  system.  Using  the  rotation  matrix  Rj^r  to  effect  the  transformation 
gives 

UJ*  =  flx  —  Clz  Sin  fly , 

Uy  —  flz  sin  flx  cos  fly  +  fly  cos  flx ,  (F10) 

uf  =  flz  cos  flx  cos  fly  —  fly  sin  flx, 

where  the  superscript  R  indicates  that  the  components  are  given  with  respect  to  the 
second  coordinate  system.  This  recovers  the  formulae  given  in  [4]  (page  609). 

For  completeness  it  is  necessary  to  comment  on  the  case  where  fly  =  ±7r/2  as  in  this  case 
there  is,  given  a  rotation  matrix  and  its  time  derivative,  indeterminancy  in  flx,  f lz,  fix 
and  flz.  If  fly  —  7t/2  then  the  angular  velocity  is 

(—fly sin  flz  \ 

fly  cos  flz  I  ,  (Fl  l) 

flz -tlx  J 


and  if  fly  =  —tt/2  the  angular  velocity  is 


—fly  sin  fl2 


u>  ~  I  fly  cos  f lz 

\  flz  +  fix 


Using  (E10)  and  (E15)  we  get  for  the  case  fly  =  7t/2 


v=(-Rn  )  ,  (F13) 

V  R'l'ij  R21 ) 

provided  R21  7^  0.  If  R21  —  0  the  last  component  of  u>  is  replaced  by  —R21/R22  which  is 
guaranteed  non-zero  in  this  case. 

Similarly,  using  (E10)  and  (E16),  we  get  for  the  case  fly  =  —  tt/2 

(  -An  \ 

w  =  .  Rn  ,  (F14) 

V  R22/R21 J 


provided  R21  7^  0.  If  R21  =  0  the  last  component  of  u>  is  replaced  by  —R21/R22  which  is 


guaranteed  non-zero  in  this  case. 


Hence  we  see  that  although  there  is  indeterminancy  in  two  of  the  Euler  angles  and  their 
time  derivatives  in  this  case,  the  components  of  the  angular  velocity  are  unique,  as  they 
should  be  given  they  reflect  a  rotational  motion  assumed  continuous. 
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